##########################################################################################
## https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
library(Scissor)
library(data.table)
library(optparse)
library(ggplot2)

##########################################################################################

option_list <- list(
    make_option(c("--single_cell_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    single_cell_file <- paste0(work_dir,"/public_ref/singleCell/pbmc_MT20_nor_PCA_50_RE0.3.Rdata")
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/singleCell"

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

single_cell_file <- opt$single_cell_file
out_path <- opt$out_pat

dir.create(out_path)

##########################################################################################
## 目前已经注释好了哪些为上皮细胞，进行重新聚类
sc_dataset <- load(single_cell_file, verbose = F)
sc_dataset_all <- pbmc_MT20_nor_PCA_50_RE0.3
##Idents函数定义要取的类别
Idents(sc_dataset_all) <- "orig.ident"   

##########################################################################################
########################################## 上皮细胞重聚类
## 提取上皮细胞
Cells.sub <- subset(sc_dataset_all@meta.data, celltype=="Epithelium")
scRNAsub <- subset(sc_dataset_all, cells=row.names(Cells.sub))
epi <- scRNAsub

## 对数据进行标准化
## 因为每个细胞测序深度不一样, 测序总reads数量不一样,所以需要进行归一化来去除测序深带来的差异,使样本间具有可比性
epi_nor <- NormalizeData(object = epi, normalization.method = 'LogNormalize',scale.factor = 10000)

## 提取在细胞间变异系数较大的基因,默认提取前2000个基因
## 鉴定表达高变基因(2000个）,用于下游分析,如PCA；
epi_nor <- FindVariableFeatures(object = epi_nor,selection.method='vst',mean.function = ExpMean,dispersion.function = LogVMR,mean.cutoff=c(0.125,3),dispersion.cutoff=c(0.5,Inf))

## 输出特征方差图
## 提取表达量变化最高的10个基因
top10 <- head(x = VariableFeatures(object=epi_nor),10)
plot1 <- VariableFeaturePlot(object = epi_nor)
plot2 <- LabelPoints(plot = plot1,points = top10,repel =T)
p <- CombinePlots(plots=list(plot1,plot2))
out_name <- paste0( out_path , "/featureVar_MT50.png"  )
ggsave(out_name, p, width = 10, height = 6)

## PCA分析
epi_nor_PCA=ScaleData(epi_nor,vars.to.regress = c("percent.mt","nCount_RNA"))
epi_nor_PCA=RunPCA(object=epi_nor_PCA,npcs=100,pc.genes=VariableFeatures(object=epi_nor_PCA))
p <- ElbowPlot(object=epi_nor_PCA,ndims=100)
out_name <- paste0( out_path , "/PCA_MT50_100.png"  )
ggsave(out_name, p, width = 10, height = 6)

###每个PC的p值分布和均匀分布
epi_nor_PCA <- JackStraw(object=epi_nor_PCA,num.replicate=100,dims=50)
epi_nor_PCA_50 <- ScoreJackStraw(object=epi_nor_PCA,dims=1:50)
p<-JackStrawPlot(object=epi_nor_PCA_50,dim=1:50)
out_name <- paste0( out_path , "/pcaJackStraw_MT50_PCA50.png"  )
ggsave(out_name, p, width = 8, height = 6)

pcSelect=50
epi_nor_PCA_50 <- FindNeighbors(object=epi_nor_PCA,dims=1:pcSelect)#计算邻接距离
## resolution设置为2提高分辨率
epi_nor_PCA_50_RE0.5 <- FindClusters(object=epi_nor_PCA_50,resolution=1)#对细胞分组，优化标准模块化

#tSNE
epi_nor_PCA_50_RE0.5<-RunTSNE(object=epi_nor_PCA_50_RE0.5,dims=1:50)#TSNE聚类

p <- TSNEPlot(object=epi_nor_PCA_50_RE0.5,label.size=4,pt.size=0.3,label=T)
out_name <- paste0( out_path , "/epi_PCA50_RE0.5_cluster.png"  )
ggsave(out_name, p, width = 6.5, height = 5.5)

p <- TSNEPlot(object=epi_nor_PCA_50_RE0.5,pt.size=0.3,group.by='orig.ident')
out_name <- paste0( out_path , "/epi_PCA50_RE0.5_sample.png"  )
ggsave(out_name, p, width = 6.5, height = 5.5)

p <- TSNEPlot(object=epi_nor_PCA_50_RE0.5,pt.size=0.3,label=T,label.size=3,group.by='celltype')
out_name <- paste0( out_path , "/epi_PCA50_RE0.5_celltype.png"  )
ggsave(out_name, p, width = 6.5, height = 5.5)

##########################################################################################

p <- FeaturePlot(object=epi_nor_PCA_50_RE0.5,features=c("MUC5AC","ATP4A","MUC6","OLFM4","CEACAM5","CEACAM6","MUC2","FABP1","CHGA","PGA3","PGA4"),cols=c('grey','red'),label=F,pt.size=0.3)
out_name <- paste0( out_path , "/featureplot2_epi_PCA50_RE0.5.png"  )
ggsave(out_name, p, width = 24, height = 18)

logFCfilter=0.5
adjPvalFilter=0.05
ALL.markers<-FindAllMarkers(object=epi_nor_PCA_50_RE0.5,
                            only.pos=F,
                            min.pct=0.25,
                            logfc.threshold=logFCfilter)                        
sig.markers <- subset(ALL.markers,(as.numeric(as.vector(ALL.markers$avg_log2FC))>logFCfilter)&(as.numeric(as.vector(ALL.markers$p_val_adj))<adjPvalFilter))
out_name <- paste0( out_path , "/epimarker_PCA50_RE0.5.xls"  )
write.table(sig.markers,file=out_name,sep='\t',row.names=F,quote=F)

## https://cloud.tencent.com/developer/article/1983799
epi_nor_PCA_50_RE0.5@meta.data$celltype[epi_nor_PCA_50_RE0.5@meta.data$seurat_clusters %in% c("4")]<-"EGC"
## MUC6
epi_nor_PCA_50_RE0.5@meta.data$celltype[epi_nor_PCA_50_RE0.5@meta.data$seurat_clusters %in% c("14" , "6" , "21")]<-"Neck"
## MUC5AC
epi_nor_PCA_50_RE0.5@meta.data$celltype[epi_nor_PCA_50_RE0.5@meta.data$seurat_clusters %in% c("3" , "8" , "5" , "9" , "2" , "0" , "16" , "23" , "1")]<-"Pit"
## FABP1
epi_nor_PCA_50_RE0.5@meta.data$celltype[epi_nor_PCA_50_RE0.5@meta.data$seurat_clusters %in% c("7" , "13" , "15" , "24" , "12")]<-"Enterocytes"
## MUC2
epi_nor_PCA_50_RE0.5@meta.data$celltype[epi_nor_PCA_50_RE0.5@meta.data$seurat_clusters %in% c("18")]<-"Goblet"
## CHGA
epi_nor_PCA_50_RE0.5@meta.data$celltype[epi_nor_PCA_50_RE0.5@meta.data$seurat_clusters %in% c("17" , "10" , "19" , "20" , "11")]<-"Endocrine"
## PGA3
epi_nor_PCA_50_RE0.5@meta.data$celltype[epi_nor_PCA_50_RE0.5@meta.data$seurat_clusters %in% c("22")]<-"Chief"

#epi_nor_PCA_50_RE0.5@meta.data$celltype[epi_nor_PCA_50_RE0.5@meta.data$seurat_clusters %in% c("0","6","7","9")]<-"ISC"

p <- TSNEPlot(object=epi_nor_PCA_50_RE0.5,pt.size=0.3,label=T,label.size=3,group.by='celltype')
out_name <- paste0( out_path , "/epi_PCA50_RE0.5_subcelltype.png"  )
ggsave(out_name, p, width = 6.5, height = 5.5)

## 输出注释好的细胞类型
out_name <- paste0( out_path , "/epi_nor_PCA_50_RE0.5.Rdata"  )
save(epi_nor_PCA_50_RE0.5,file=out_name)
